
struct Pos
	{
	int x, y;
	};


struct Z
	{
	unsigned z1, z2, z3, z4;
	};

// A and C are constants  
//__device__ unsigned lcgStep(unsigned &z, unsigned A, unsigned C)  
//	{  
//	return z = (A * z + C);  
//	}

#define lcgStep(z, A, C)   (z = (A * z + C))  

// S1, S2, S3, and M are all constants, and z is part of the  
// private per-thread generator state.  
__device__ unsigned tausStep(unsigned &z, int S1, int S2, int S3, unsigned M)  
	{  
//	const unsigned b=(((z << S1) ^ z) >> S2);  
	return z = (((z & M) << S3) ^ (((z << S1) ^ z) >> S2));  
	}

__device__ unsigned hybridTausU(Z & z)  
	{  
	// combined period is lcm(p1,p2,p3,p4)~ 2^121  
	return   
		tausStep(z.z1, 13, 19, 12, 4294967294ul) ^  // p1=2^31-1  
		tausStep(z.z2, 2, 25, 4, 4294967288ul) ^    // p2=2^30-1  
		tausStep(z.z3, 3, 11, 17, 4294967280ul) ^   // p3=2^28-1  
		lcgStep(z.z4, 1664525, 1013904223ul)        // p4=2^32  
		;  
	}  

__device__ unsigned sq_dist(const Pos & p1, const Pos & p2)
	{
	return (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);
	}


/** Convert a 2-dimensional coordinate into a z-value. */
__device__ inline unsigned point2z(const int x, const int y)
	{
	unsigned z = 0;
	
	// coords are half the size of z-Values
	for(unsigned j=0; j<8*sizeof(unsigned)/2; j++)
		{
		const int mask = 1 << j;
		
		z |= (unsigned(x & mask) << j) | 
			 (unsigned(y & mask) << (j+1));
		}
	
	return z;
	}

__device__ inline void z2point(unsigned z, int & x, int & y)
	{
	x = 0; y = 0;

	unsigned mask = 1;
	
	for (unsigned j=0; j<8*sizeof(unsigned)/2; j++)
		{
		x |= unsigned(z & mask);
		z >>= 1;
		y |= unsigned(z & mask);
		mask <<= 1;
//			point.x |= coord_type((z & key_type(1)<<j*2)   >> j);
//			point.y |= coord_type((z & key_type(1)<<j*2+1) >> j+1);
		}
	}

__device__ inline int sq_dist(int x1, int y1, int x2, int y2)
	{
	return (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2);
	}

struct SortData
	{
	unsigned z;
	int i;
	};

template<int CHUNK_SZ, SortData * SOURCE, SortData * TARGET>
__device__ void merge()
	{
	const int chunk = threadIdx.x / CHUNK_SZ;
	const int odd = chunk % 2;
	const int other = chunk + (1 - 2*odd);

	const int rank = threadIdx.x % CHUNK_SZ;
	const int self = SOURCE[chunk*CHUNK_SZ + rank].z + odd;

	int o_rank = rank;
	while (o_rank < (CHUNK_SZ) && SOURCE[other*CHUNK_SZ+o_rank].z < self)
		o_rank++;

	while(o_rank > 0 && SOURCE[other*CHUNK_SZ+o_rank-1].z >= self)
		o_rank--;

	TARGET[(chunk/2) * 2 * CHUNK_SZ + rank + o_rank].z = self - odd;
	}

//__shared__ unsigned s_coll[256];
__shared__ SortData s_sort[256];
__shared__ SortData s_sort_tmp[256];
__shared__ int sort_idx[256];
__shared__ Pos s_pos[256];


__global__ void scalarProdGPU(Pos * iPos, Z * rngz, unsigned * coll)
	{
	const int i = blockDim.x * blockIdx.x + threadIdx.x;
	const int t = threadIdx.x;

//	Pos pos = iPos[i];
	Z z = rngz[i];
	
	// no safety against mismatch of data size <-> thread number
	Pos pos = iPos[i];
	unsigned tcoll = coll[i];

	pos.x = (pos.x + hybridTausU(z) % 5 - 2 + 1024) % 1024;
	pos.y = (pos.y + hybridTausU(z) % 5 - 2 + 1024) % 1024;

	const unsigned zValue = s_sort[t].z = point2z(pos.x, pos.y);
	s_sort[t].i = t;
	
	//s_pos[t] = pos;
	__syncthreads();
	
	merge<1, s_sort, s_sort_tmp>();
	__syncthreads();
	merge<2, s_sort_tmp, s_sort>();
	__syncthreads();
	merge<4, s_sort, s_sort_tmp>();
	__syncthreads();
	merge<8, s_sort_tmp, s_sort>();
	__syncthreads();
	merge<16, s_sort, s_sort_tmp>();
	__syncthreads();
	merge<32, s_sort_tmp, s_sort>();
	__syncthreads();
	merge<64, s_sort, s_sort_tmp>();
	__syncthreads();
	merge<128, s_sort_tmp, s_sort>();
	__syncthreads();

	// where position of originally ith threat ended up
//	sort_idx[s_sort[t].i] = t;
	z2point(s_sort[t].z, s_pos[t].x, s_pos[t].y);

	__syncthreads();

	const int visibility = 10;
	const int sq_vis = visibility * visibility;

	int lx = max(0, s_pos[t].x - visibility);
	int hx = min(1023, s_pos[t].x + visibility);

	int ly = max(0, s_pos[t].y - visibility);
	int hy = min(1023, s_pos[t].y + visibility);

	const unsigned lz = point2z(lx, ly);
	const unsigned hz = point2z(hx, hy);

	int lidx = t - 1;
	while(lidx >= 0 && s_sort[lidx].z >= lz)
		{
		if (sq_dist(s_pos[t].x, s_pos[t].y, s_pos[lidx].x, s_pos[lidx].y) < 
			sq_vis)
			tcoll ++;

		lidx--;
		}
	
	int idx = t + 1;
	while(idx<256 && s_sort[idx].z <= hz)
		{
		if (sq_dist(s_pos[t].x, s_pos[t].y, s_pos[idx].x, s_pos[idx].y) < 
			sq_vis)
			tcoll ++;

		idx++;
		}

// O(n^2):

//	for (int x=0; x<256; x++)
//		if (sq_dist(pos.x, pos.y, s_pos[x].x, s_pos[x].y) < 100)
//			s_coll[t]++;	

	coll[i] = tcoll;
	iPos[i] = pos;
	rngz[i] = z;
	}

